11. Estimadores adaptivos parte I¶

Hasta ahora hemos estudiando sistemas lineales donde:

sus coeficientes quedan fijos luego del diseño y son constantes en el tiempo
hacen supuestos sobre los estadísticos de la señal/ruido

¿Qué podemos hacer si

no podemos hacer supuestos sobre los estadísticos?
los estadísticos de la señal/ruido cambian en el tiempo?
estamos en un escenario donde los datos llegan continuamente (data streaming)?

Cuando los estadísticos cambian en el tiempo decimos que la señal es no estacionaria.

En estos casos necesitamos un estimador de tipo adaptivo, es decir sistemas y filtros cuyos coeficientes se pueden adaptar a medida que llegan nuevos datos. Estos estimadores se diseñan de acuerdo a un método de optimización que es online

La siguiente figura muestra algunos ejemplos de aplicaciones de sistemas adaptivos

El método de optimización online que utilizaremos principalmente en este curso es el gradiente descendente estocástico. Revisemos a continuación los fundamentos.

11.1. Gradiente descendente¶

Sea un vector de pesos \(w\) de largo \(L+1\) que guarda los coeficientes de un estimador

Sea ahora una función de costo que mapea el vector de pesos a un número real

\[ J(w): \mathbb{R}^{L+1} \to \mathbb{R} \]

La función de costo debe ser tal que a menor \(J\) menor sea el error del estimador

Para entrenar un estimador o filtro adaptivo se tienen los siguientes pasos conceptuales

Partimos de una solución inicial \(w_0\)
Modificamos iterativamente \(w\) tal que \(J(w_{t+1}) < J(w_t)\)
Nos detenemos al cumplir un cierto criterio

Para modificar iterativamete y eficientemente los pesos utilizaremos la regla del gradiente descendente (GD)

\[ w_{t+1} = w_t - \mu \frac{dJ(w)}{dw}, \]

donde \(\mu\) se conoce como tasa de aprendizaje o “paso”

Imaginemos \(J\) como una superficie de \(L+1\) dimensiones
En cada punto el gradiente negativo de \(J\) nos indica hacia donde está el descenso más abrupto
La tasa \(\mu\) nos da el largo del salto entre \(w_t\) y \(w_{t+1}\)

Observando la expansión de Taylor de primer orden de \(J\) en \(w_{t}\)

\[\begin{split} \begin{align} J(w_{t+1}) &= J(w_t) + \frac{dJ(w_t)}{dw} (w_{t+1} - w_{t}) \nonumber \\ &= J(w_t) -\mu \left \| \frac{dJ(w_t)}{dw} \right \|^2 \leq J(w_t) \nonumber \end{align} \end{split}\]

es decir que usando la regla GD con \(\mu>0\) y asumiendo que \(J\) es convexo entonces se cumple que \(J\) siempre decrece monotónicamente.

La siguiente gráficas interactivas muestran una superficie de costo no convexa para un parámetro unidimensional. Cada punto representa una solución que parte desde una posición inicial distinta. Las flechas corresponden a la derivada multiplicada por la tasa de aprendizaje.

Estudie la evolución de las tres soluciones en cada caso. En primer lugar se utiliza \(\mu=0.05\)

J = lambda w : (w-1)**2 + 0.2*np.sin(2*np.pi*w) #  Función de costo
gradJ = lambda w : 2*(w-1) + 0.2*2*np.pi*np.cos(2*np.pi*w) # Gradiente
mu = 0.05 # Tasa de aprendizaje
iteraciones = 15
wt = np.zeros(shape=(iteraciones, 3))
wt[0, :] = np.array([0.05, 0.4, 1.9]) # Solución inicial
w_plot = np.linspace(0, 2, num=100)

for k in range(1, iteraciones):
    wt[k, :] = wt[k-1, :] - mu*gradJ(wt[k-1, :])

loss_surface = hv.Curve((w_plot, J(w_plot)), 'w', 'J')
hMap = hv.HoloMap(kdims='Iteración')
for k in range(iteraciones):
    dots = hv.Points((wt[k, :], J(wt[k, :]))).opts(size=10, color='k')
    mag = mu*gradJ(wt[k, :])
    angle = np.pi/2 - np.sign(-mag)*np.pi/2
    mag = np.abs(mag)
    arrows = hv.VectorField((wt[k, :], J(wt[k, :]), angle, mag)).opts(pivot='tail', 
                                                                      magnitude=hv.dim('Magnitude'), 
                                                                      rescale_lengths=False)
    hMap[k] =  dots * arrows 
    
loss_surface * hMap

Advertencia

Dependiendo de donde partimos la solución final es distinta. El gradiente descedente puede quedarse “atorado” en un mínimo local o en un punto silla

Ahora observe como evolucionan las tres soluciones con \(\mu=0.5\), es decir 10 veces más grande que el caso anterior

J = lambda w : (w-1)**2 + 0.2*np.sin(2*np.pi*w) #  Función de costo
gradJ = lambda w : 2*(w-1) + 0.2*2*np.pi*np.cos(2*np.pi*w) # Gradiente
mu = 0.5 # Tasa de aprendizaje
iteraciones = 15
wt = np.zeros(shape=(iteraciones, 3))
wt[0, :] = np.array([0.05, 0.4, 1.9]) # Solución inicial
w_plot = np.linspace(0, 2, num=100)

for k in range(1, iteraciones):
    wt[k, :] = wt[k-1, :] - mu*gradJ(wt[k-1, :])

loss_surface = hv.Curve((w_plot, J(w_plot)), 'w', 'J')
hMap = hv.HoloMap(kdims='Iteración')
for k in range(iteraciones):
    dots = hv.Points((wt[k, :], J(wt[k, :]))).opts(size=10, color='k')
    mag = mu*gradJ(wt[k, :])
    angle = np.pi/2 - np.sign(-mag)*np.pi/2
    mag = np.abs(mag)
    arrows = hv.VectorField((wt[k, :], J(wt[k, :]), angle, mag)).opts(pivot='tail', 
                                                                      magnitude=hv.dim('Magnitude'), 
                                                                      rescale_lengths=False)
    hMap[k] =  dots * arrows 
    
loss_surface * hMap

Advertencia

Si la tasa de aprendizaje es muy alta, los pasos son muy largos y podríamos no converger a un punto estacionario

Los ejemplos anteriores nos han mostrado algunas de las limitaciones del algoritmo de gradiente descendente. Es importante tenerlas en cuenta cuando lo utilicemos en nuestras aplicaciones

11.2. Gradiente descendente en el filtro de Wiener¶

Para el filtro de Wiener teníamos que

\[ J(h) = \sigma_d^2 - 2 \textbf{h}^T R_{ud} + \textbf{h}^T R_{uu} \textbf{h}, \]

por ende

\[ \frac{dJ(h)}{dh} = -2 R_{ud} + 2 R_{uu} \textbf{h} \]

y finalmente

\[ \textbf{h}_{t+1} = \textbf{h}_{t} (I - 2 \mu R_{uu}) + 2\mu R_{ud} \]

En este caso la condición para una convergencia estable es

\[ 0 < \mu < \frac{1}{\lambda_{\text{max}}}, \]

donde \(\lambda_{\text{max}}\) es valor propio más grande de \(R_{uu}\)

(La prueba de esto puede encontrarse en Haykin, “Adaptive filter theory”, Sección 4.3)

11.3. Gradiente descendente estocástico (SGD)¶

El filtro de Wiener es óptimo pero no adaptivo:

Requiere de \(N\) muestras de \(u\) y \(d\) para estimar \(R_{ud}\) y \(R_{uu}\)
Los pesos se adaptan luego de haber presentado las \(N\) muestras: Es una estrategia de tipo batch
Asume que la señal es estacionaria

Si nuestros son no estacionarios significa que debemos adaptar el filtro a medida que nuevas muestras son observadas . Para lograr esto podemos usar la versión estocástica del GD: SGD

En SGD:

los pesos se adaptan luego de haber presentado una sola muestra o un conjunto pequeño de muestras (mini-batch)
no hay garantía de llegar al óptimo en un problema convexo, pero es más eficiente computacionalmente que GD

El siguiente esquema muestra una comparación entre la trayectoria de \(w\) cuando se usa GD (negro) y SGD (rojo). En general la trayectoria de SGD será más ruidosa y también podría requerir más pasos, pero cada paso es mucho más económico

11.4. Algoritmo Least Mean Square (LMS)¶

Podemos extender el filtro de Wiener al caso no-estacionario usando SGD, el resultado es un algoritmo simple que además es robusto: El algoritmo LMS

Fue fue inventado en 1960 por Bernard Widrow y Ted Hoff
A diferencia del filtro de Wiener no se requiere conocimiento estadístico del proceso. Tampoco se requiere calcular e invertir la matriz de correlación
El algoritmo LMS se ajusta o entrena de manera recursiva y online

Consideremos la función de costo estocástica para la arquitectura FIR que utilizamos para el filtro de Wiener

\[\begin{split} \begin{align} J^s_n(\textbf{w}) &= e_n^2 \nonumber \\ &= (d_n - y_n)^2 \nonumber \\ &= (d_n - \textbf{w}^T \textbf{u}_n )^2 \nonumber \\ &= (d_n - \sum_{k=0}^{L} w_{n, k} u_{n-k} )^2 \nonumber \end{align} \end{split}\]

donde definimos \(\textbf{u}_n = [u_n, u_{n-1}, \ldots, u_{n-L}]\).

Nota

A diferencia del filtro de Wiener no aplicamos el valor esperado al error cuadrático. Se usa el error cuadrático instantaneo

Para continuar calculamos el gradiente en función del peso \(w_{n, k}\)

\[ \frac{d J^s_n (\textbf{w})}{d w_{n, k}} = - 2 e_n u_{n-k} \]

Luego, usando la regla SGD llegamos a

\[ w_{n+1, k} = w_{n, k} + 2 \mu e_n u_{n-k}, k=0, 1, \ldots, L \]

y que en forma matricial es

\[\begin{split} \begin{align} \textbf{w}_{n+1} &= \textbf{w}_{n} + 2 \mu e_n \textbf{u}_{n}\nonumber \\ &= \textbf{w}_{n} + 2 \mu (d_n - \textbf{w}_{n}^T \textbf{u}_{n}) \textbf{u}_{n}, \nonumber \end{align} \end{split}\]

que se conoce como la regla de Widrow-Hoff

Importante

El algoritmo LMS estima el error instantaneo y actualiza los pesos recursivamente

La complejidad de este algoritmo es \(L+1\).

11.4.1. Convergencia del algoritmo LMS (Haykin 6.5)¶

El algoritmo LMS tiende en la media al valor óptimo

\[ \mathbb{E}[\textbf{w}_n] \to \textbf{w}^* \]

para \(n\to \infty\)

Además convergence en la media cuadrada: La varianza de \(\textbf{w}_n - \textbf{w}^*\) tiene al valor mínimo de \(J\) para \(n\to \infty\)

Esto se cumple si

\[ 0 < \mu < \frac{2}{\text{Tr}[R_{uu}]} \]

donde \(R_{uu} = \mathbb{E}[\textbf{u}_n \textbf{u}_n^T ]\) es la matriz de autocorrelación y \(\text{Tr}[]\) el operador que calcula la traza de una matriz

11.4.2. Algoritmo Normalized LMS (NLMS)¶

Tenemos la siguiente regla iterativa

\[\begin{split} \begin{align} \textbf{w}_{n+1} &= \textbf{w}_{n} + 2 \mu (d_n - \textbf{w}_{n}^T \textbf{u}_{n}) \textbf{u}_{n} \nonumber \\ & = \textbf{w}_{n} + \Delta \textbf{w}_n \nonumber \end{align} \end{split}\]

que se puede interpretar graficamente como

(donde \(\textbf{x}(k)\) y \(\textbf{w}(k)\) corresponden a \(\textbf{u}_n\) y \(\textbf{w}_n\) en nuestra notación, respectivamente)

Nota

Los cambios en el vector de peso \(\Delta \textbf{w}_n\) son paralelos a \(\textbf{u}_{n}\). Además estos cambios podrían estar dominados por

\[ \max \textbf{u}_{n} = [u_n, u_{n-1}, \ldots, u_{n-L}] \]

El algoritmo Normalized LMS (NLMS) corrige este problema ponderando por la varianza de \(\textbf{u}_{n}\)

\[ \textbf{w}_{n+1} = \textbf{w}_{n} + 2 \mu (d_n - \textbf{w}_{n}^T \textbf{u}_{n}) \frac{\textbf{u}_{n}}{\left(\|\textbf{u}_{n}\|^2 + \delta\right)} \]

donde la constante \(\delta\) es un valor pequeño que se usa para evitar divisiones por cero. En lo que sigue usaremos NLMS para revisar algunas aplicaciones

11.5. Implementación del filtro NLMS en Python¶

Podemos implementar las ecuaciones del filtro NLMS como se muestra a continuación

class Filtro_NLMS:
    
    def __init__(self, L, mu, delta=1e-6, winit=None):
        self.L = L
        self.w = np.zeros(shape=(L+1, ))
        self.mu = mu
        self.delta = delta
        
    def update(self, un, dn):
        # Asumiendo que un = [u[n], u[n-1], ..., u[n-L]]
        unorm = np.dot(un, un) + self.delta
        yn = np.dot(self.w, un)
        self.w += 2*self.mu*(dn - yn)*(un/unorm)
        return yn

El filtro recibe como entrada el orden \(L\) y la tasa de aprendizaje \(\mu\)
Se asume un vector cero para los pesos iniciales, pero también en la práctica podemos partir de una solución anterior si esta existiera
Para actualizar el vector de pesos es necesario entregar el vector \(\textbf{u}_n \in \mathbb{R}^{L+1}\) y la salida deseada \(d_n \in \mathbb{R}\). La función update retorna la salida predicha por el filtro \(y_n = w_n^T \textbf{u}_n \)

A continuación probaremos este filtro con una aplicación conocida como Adaptive line enhancement (ALE). ALE se refiere a un sistema adaptivo para eliminar ruido blanco aditivo de una señal. El sistema aprende un filtro pasabanda en torno a la frecuencia de interés

En ALE usamos como señal deseada

\[ d_n = u_n = \textbf{u}_n[0] \]

El valor predicho por el filtro será la señal \(u\) pero libre de ruido blanco. Esto se debe a que el ruido blanco no tiene correlación y por ende el filtro adaptivo no lo puede predecir

# Digamos que u = s + n
# El objetivo es limpiar u para obtener s
# s es una señal determínista y n es ruido blanco

Fs, f0 =  100, 5
t = np.arange(0, 4, 1/Fs)
s = np.sin(2.0*np.pi*t*f0)
n = 0.5*np.random.randn(len(t)) 
s[t>2.0] += 5  # Simulemos un cambio abrupto en la media de la señal
#s += s*(0.5 + 0.5*np.cos(2.0*np.pi*t/2))  # Tremolo (AM)
u = s + n

A diferencia de un filtro estático (como el filtro de Wiener) es posible filtrar incluso ante cambios bruscos en la señal.

Estudie como cambia el resultado del filtro con distintos valores de \(\mu\)

L = 20
u_preds = {}

for mu in np.logspace(-2, 0, num=10):
    myfilter = Filtro_NLMS(L=L, mu=mu)
    u_preds[mu] = np.zeros(shape=(len(u),))
    for k in range(L+1, len(u)):
        u_preds[mu][k] = myfilter.update(u[k-L-1:k][::-1], u[k])

hMap = hv.HoloMap(kdims='mu')
for mu, u_pred in u_preds.items():
    s1 = hv.Curve((t, s), 'Tiempo', 'Señal', label='Limpia')
    s2 = hv.Scatter((t, u), 'Tiempo', 'Señal', label='Contaminada')
    s3 = hv.Curve((t, u_pred), 'Tiempo', 'Señal', label='Filtrada')
    hMap[mu] = hv.Overlay([s1, s2, s3]).opts(hv.opts.Overlay(legend_position='top'), 
                                             hv.opts.Curve(ylim=(-5, 10), height=350))
hMap

Importante

La tasa de aprendizaje \(\mu\) controla la velocidad de adaptación. Pero una tasa demasiado grande provoca que el filtro sea inestable. En general el valor óptimo de \(\mu\) depende del problema y del valor de \(L\)

La siguiente figura muestra la respuesta en frecuencia del filtro en función del tiempo para \(\mu=0.02\)

Observe como a medida que se adapta el filtro se concentra en la frecuencia fundamental de la señal, que en este caso es 5 Hz

L = 20
u_preds = {}
myfilter = Filtro_NLMS(L=L, mu=0.02)
H_history = np.zeros(shape=(512, len(u))) 
for k in range(L+1, len(u)):
    myfilter.update(u[k-L-1:k][::-1], u[k])
    fk, Hk = scipy.signal.freqz(b=myfilter.w, a=1, fs=Fs)
    H_history[:, k] = np.abs(Hk)

hv.Image((t, fk, H_history), kdims=['Tiempo [s]', 'Frecuencia [Hz]']).opts(cmap='Blues') 

11.6. Comparación entre Filtro de Wiener/GD y algoritmo LMS/SGD¶

Supuestos: Wiener requiere un ambiente estacionario lo cual nos permite calcular \(R_{uu}\) y \(R_{ud}\). En LMS la señal puede ser no estacionaria.
Aprendizaje: En el filtro de Wiener el aprendizaje es determinista. En LMS el aprendizaje viene promediando a nivel de los estimadores de \(w\). En LMS el aprendizaje es estadístico.
Optimalidad: Wiener es óptimo en cambio LMS es sub-óptimo (localmente óptimo). LMS tiende a la solución de Wiener
Costo: LMS se actualiza online y tiene costo \(L\). Wiener se entrena offline y tiene costo \(L^2\)

A continuación se muestra un diagrama que compara el filtro de Wiener y el algoritmo LMS

INFO183 Análisis de sistemas lineales